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A Riemann problem with prescribed initial conditions will produce one of three possible wave 
patterns corresponding to the propagation of the different discontinuities that will be produced 
once the system is allowed to relax. In general, when solving the Riemann problem numerically, 
the determination of the specific wave pattern produced is obtained through some initial guess 
which can be successively discarded or improved. We here discuss a new procedure, suitable for 
implementation in an exact Riemann solver in one dimension, which removes the initial ambi- 
guity in the wave pattern. In particular we focus our attention on the relativistic velocity jump 
between the two initial states and use this to determine, through some analytic conditions, the 
wave pattern produced by the decay of the initial discontinuity. The exact Riemann problem is 
then solved by means of calculating the root of a nonlinear equation. Interestingly, in the case 
of two rarefaction waves, this root can even be found analytically. Our procedure is straightfor- 
ward to implement numerically and improves the efficiency of numerical codes based on exact 
Riemann solvers. 



1. Introduction 

The general Riemann problem is based on the calculation of the one-dimensional temporal 
evolution of a fluid which, at some given initial time, has two adjacent states characterized by 
different values of uniform velocity, pressure and density. From a mathematical point of view, 
the Riemann problem is an initial value problem for a hyperbolic system of partial differential 
equations, with initial conditions characterized by a discontinuity between the two initial states. 
These initial conditions establish the way in which the discontinuity will decay after removal of 
the barrier separating the initial "left" and "right" states. The schematic evolution of a general 
Riemann problem can be represented as (Marti & Miiller 1994) 

LW^LXR*W^R, (1.1) 

where W denotes a shock or a rarefaction wave that propagates towards the left («— ) or the right 
{-^) with respect to the initial discontinuity, L and R are the known initial left and right states, 
while and i?* are the new hydrodynamic states that form behind the two waves propagating 
in opposite directions. These waves are separated by a contact discontinuity C and therefore have 
the same values of the pressure and velocity, but different values of the density|f[ The Riemann 
problem is said to be solved when the velocity, pressure and density in the new states and 
i?* have been computed, as well as the positions of the waves separating the four states. The 
solution of the one dimensional Riemann problem in relativistic hydrodynamics was discussed 

f Note that the contact discontinuity might be also trivial so that the density may not be discontinuous. 
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in the general case by Marti & Miiller (1994) and the reader is referred to their work for further 
details (see also Pons et al. 2000, for the extension to multidimensions). 

The numerical solution of a Riemann problem is the fundamental building block of hydro- 
dynamical codes based on Godunov-type finite difference methods. In such methods, the com- 
putational domain is discretized and each interface between two adjacent grid-zones is used to 
construct the initial left and right states of a "local Riemann problem". The evolution of the hy- 
drodynamical equations is then obtained through the solution across the computational grid of 
the sequence of local Riemann problems set up at the interfaces between successive grid-zones 
(see Marti & Miiller 1996, but also Godunov 1959 and Colella & Woodward 1984). The "core" of 
each of these local Riemann problems consists of determining the fluid pressure across the con- 
tact discontinuity, which can be calculated by imposing the continuity of the normal component 
of the fluid velocity across C 

vlAp*) = '"rAp*) ■ (1-2) 



In general, (1.2) is a nonlinear algebraic equation in the unknown pressure and requires a 
numerical solution. Depending on the different wave patterns forming after the decay of the 
discontinuity, a different nonlinear equation will need to be solved|[ This initial "ambiguity" in 
the wave pattern produced corresponds to the fact that the interval in pressure bracketing the 
solution is not known a priori. In practice this lack of information is compensated by the use 
of efficient numerical algorithms which, via a process of trial and error, determine the correct 
wave pattern and then proceed to the solution of the corresponding nonlinear equation (Marti, & 
MuUer 1999). 

In this paper, we show that the relativistic expression for the relative velocity between the 
two initial states is a function of the unknown pressure and so a new procedure for numeri- 
cally solving the exact Riemann problem can be proposed in which the pressure is no longer 



obtained by the solution of equation (1.2). Rather, p» is calculated by equating the relativistic 
invariant expression for the relative velocity between the two initial states with the value given 
by the initial conditions. 

When compared to equivalent approaches, our exact Riemann solver has some advantages. 
Firstly, we can remove the ambiguity mentioned above and determine the generated flow pattern 
by simply comparing the relative velocity between the two initial states with reference values 
built from the initial conditions of the Riemann problem. Doing so provides immediate infor- 
mation about which of the nonlinear equations (one for every wave pattern) needs to be solved. 
Secondly, by knowing the wave pattern we can produce an immediate bracketing of the solu- 
tion. Doing so gives improved efficiency in the numerical root finding procedure. Finally, for one 
of the wave patterns (i.e. for two rarefaction waves moving in opposite directions) our method 
provides the solution of the relativistic Riemann problem in a closed analytic form. 

Because of its simplicity, the numerical implementation of our method is straightforward and 
can be accomplished with a much smaller number of lines of code. When compared with other 
exact Riemann solvers (e. g. Marti & Miiller 1999) it has also proved to be computationally more 
efficient. In particular, when solving a generic hydrodynamical problem (in which one solves for 
very simple Riemann problems) the approach proposed here brackets the solution very closely 
and this produces substantial computational improvements of up to 30%. For the cases discussed 
here however, where very strong shocks are considered, the speed-up is below 10%. 

This paper briefly introduces our idea and is organized as follows: in Section ^ we define the 
mathematical set-up for the formulation of the Riemann problem in relativistic hydrodynamics. 

f This approach is usually referred to as an "exact" Riemann solver to distinguish it from the family 
of so called "approximate" Riemann solvers, where the system of equations to be solved is reduced to 
quasi-linear form, thus avoiding any iterative procedure. See, for example, the approximate Riemann solver 
of Roe (1981). 
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In Section || we determine the basic relativistic expressions linking the velocities ahead of and 
behind a shock or a rarefaction wave. These expressions will be used repeatedly in the following 
Sections where will write explicit criteria for the occurrence of the three possible wave 
patterns. Section ^ discusses how the criteria can be used for making an efficient numerical 
implementation of a Riemann solver and briefly presents a comparison of performances with 
more traditional algorithms. Finally, Section |] presents our conclusions. Throughout, we use a 
system of units in which c = 1. 



2. Setting up the Relativistic Riemann Problem 

Consider a perfect fluid with four-velocity ~ Wil^ v, 0, 0) with W = {I — w^)^^/^ being 
the Lorentz factor and with a stress-energy tensor 

T"" ^ (e + p)u''u'' + pT]""" = phu^u" + pTj"'' , (2.1) 

where ^ = 0, . . . , 3, 77'"' ~ diag(— 1, 1,1,1) and e, p are the proper energy density and pressure, 
respectively. Assume moreover that the fluid obeys a polytropic equation of state 

p = fc(s)p'r = (7 - l)pe , (2.2) 

where p is the proper rest mass density, 7 is the adiabatic index, and k{s) is the polytropic 
constant, dependent only on the specific entropy s (the latter is generally assumed to be different 
in the two initial states). Straightforward expressions can be written relating e and p to the specific 
enthalpy h and to the specific internal energy e of the fluid 

e^p{l + e) = p+^, (2.3) 
h=l + e + P=l + P(^] . (2.4) 



P P \-f 

Consider now the fluid to consist of an initial left state (indicated with an index 1) and an initial 
right state (indicated with an index 2), each having prescribed and different values of uniform 
pressure, density and velocity|t[ The discontinuity between the two states which has been con- 
structed in this way will then decay producing one of the three wave patterns listed below 

( i) two shock waves, one moving towards the initial left state, and the other towards the initial 
right state: LS^L^,CR^,S^R. 

( ii) one shock wave and one rarefaction wave, the shock moving towards the initial right state, 
and the rarefaction towards the initial left state (or viceversa if p2 > pi): LTZ<^L.XR*S^R. 

( Hi) two rarefaction waves, one moving towards the initial left state, and the other towards the 
initial right state: LTZ^L^^CR^TZ^R. A special case of this wave pattern is produced when the 
two rarefaction waves leave a vacuum region behind them. 

The basis of our approach lies in the possibility of determining "a priori" which of these three 
wave patterns will actually result by simply comparing the value of the special relativistic relative 
velocity between the initial left and right states 

V\2 = , (2.5) 

with the values of the limiting relative velocities for the occurrence of the wave patterns (i)-(iii). 

t Note that hereafter we will consider the fluid to be of the same type in the two initial states. In prin- 
ciple, in fact, the fluid in the two initial states might be governed by two different equations of state or by 
polytropic equations of state with different polytropic indices. The formulation of the problem in that case 
is equivalent to the one presented here for a single type of fluid but particular attention must be paid to 
distinguishing the different components in the relevant expressions. 
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In this respect, our approach represents the relativistic generalization of a similar analysis pro- 
posed in Newtonian hydrodynamics by Landau and Lifshitz (1987) (see also Gheller 1997 for a 
numerical implementation). Mathematically, the occurrence of the different cases can be distin- 
guished because the relative velocity|| ( |2.5| ) is a continuous monotonic function of the pressure 

across the contact discontinuity. This is shown in Fig. |l] where we have plotted vi2 as a func- 
tion of and it is proved mathematically in Appendix A. Note that the curve shown in Fig. [l| is 
effectively produced by the continuous joining of three different curves describing the relativistic 
relative velocity for the three different wave patterns corresponding respectively to two shocks 
(2S, indicated as a dashed line), one shock and one rarefaction wave (SR, indicated as a dotted 
line), and two rarefaction waves (2R, indicated as a continuous line). The joining of the differ- 
ent curves occurs at precise values of the relative velocity which we have indicated as (z;i2)2S' 
{vi2)sj^, {vi2)2R ™d which depend uniquely on the initial conditions of the two unperturbed 
states. These three values of the relative velocity mark the extremes of three different intervals 
on the vertical axis and correspond to the three different cases (i)-(iii). 

Within this framework then, it is sufficient to compare the relative velocity between the initial 
states (2^)withthe calculated values of the limiting relative velocities ('512)23' (''^12)51?' ('^i2)2fi' 
to determine, in advance, which wave pattern will be produced. This, in turn, determines the 
correct equation to solve and the correct bracketing of the solution. In the following Sections we 
will discuss in detail how to derive analytic expressions for the limiting relative velocities and 
use them efficiently within a numerical code. 



3. Limiting Relative Velocities: Fundamental Expressions 



The expression for the relative velocity ( 2.5 ) is a relativistic invariant, but its calculation can be 
considerably simplified when performed in an appropriate reference frame. In practice, we will 
consider each of the three wave patterns ( i)-( Hi) as being composed of two "discontinuity fronts" 
(two shocks, two rarefaction waves, or one of each) moving in opposite directions and separated 
by a region where a contact discontinuity is present. In the case of a shock, in particular, it is 
useful to use a reference frame comoving with the shock front, in which the relativistic expression 
for the relative velocities ahead of (a) and behind (b) the shock takes the form (Taub 1978) 



Vab 



Vb 



1 - VaVb 



{Pb - Pa){eb - eg) 
{Ca + Pb){eb + Pa) 



(3.1) 



In the case of a rarefaction wave, on the other hand, it is more convenient to use the Eulerian 
frame in which the initial states are measured and which has one of the axes aligned with the 
direction of propagation of the wave front. In such a frame, the flow velocity at the back of a 
rarefaction wave (i.e. behind the tail of the rarefaction wave) can be expressed as a function of 
pressure at the back of the wave as 



Vb = 



{l+Va)A±{pb)-il-Va) 
{I + Va)A±{pb) + {1 - Va) 



The quantity A±{p) in (3.2) is defined as (Marti & Muller 1994) 



A±{p) 



(7-l)i/2-c,(p) 



(7-l)i/2 + c,(p) 



±2/(7-1)^/" 



(3.2) 



(3.3) 



with the ± signs corresponding to rarefaction waves propagating to the left ( TZ^ ) and to the 
right (TZ^) of the contact discontinuity, respectively. The quantity Cs{p) in (3.3) is the local 



J For compactness we will hereafter refer to as "relative velocity" the relativistic invariant expression. 
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Figure 1 . Relative velocity between the two initial states 1 and 2 as a function of the pressure at the 
contact discontinuity. Note that the curve shown is given by the continuous joining of three different curves 
describing the relative velocity corresponding respectively to two shocks (a dashed line), one shock and 
one rarefaction wave (dotted line), and two rarefaction waves (continuous line). The joining of the curves 
is indicated with filled dots. The small inset on the right shows a magnification for a smaller range of p» 
and we have indicated with filled squares the limiting values for the relative velocities (^12)23' (^12)3^' 
("12)27; 



sound speed which, for a polytropic equation of state, can be written as 



' 7(7 

(7 - 1)P + IP 



(3.4) 



We can now use expression (3.2) to write an invariant expression for the relative velocity 
across a rarefaction wave (i.e. the relative velocity of the fluid ahead of the rarefaction wave and 
behind the tail of the rarefaction wave) as 



Vab = 



1 + A±{pb) 



(3.5) 



Expressions (3.1) and (3.5) are not yet in a useful form since they cannot be combined to give 



(2.5). However, we can exploit the fact they are relativistic invariants to evaluate them in the rest 
frame of the contact discontinuity. The latter is, by definition, comoving with the fluid behind 
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the discontinuity fronts and within such a frame we can set Vb = and use equations (3.1), (3.5) 
to obtain an explicit expression for the velocity ahead of the wave front. The expressions for the 
velocities ahead of both discontinuity fronts obtained in this way can then be combined so as to 
evaluate expression (2.5). 

In the following Sections we will apply the procedure outlined above to derive the relative 
velocity between the left and right states for the three different wave patterns (i)-(iii). We will 
denote with indices 3 and 3' the quantities evaluated in the left (L,) and right (i?*) states behind 
the discontinuity fronts, and where ~ py, vj, — vy, and ^ py. Note also that in all of the 
different cases we will assume that pi > P2, and take the positive a;-direction as being the one 
from the region 1 to the region 2 (An alternative opposite choice is also possible.). 



4. LS^LXR*S^R: Two Shock Fronts 

We start by considering the wave pattern produced by two shocks propagating in opposite 
directions (see Fig. This situation is characterized by a value of the pressure downstream of the 
shocks which is larger than the pressures in the unperturbed states, i.e. p^ > pi > P2- Applying 



equation (3.1 ) to the shock front moving towards the left and evaluating it in the reference frame 



of the contact discontinuity, we can write the velocity ahead of the left propagating shock as 



Vl 



(P3 -Pi)(e3 - ei) 
(ei +P3)(e3 +Pi) 



(4.1) 



Similarly, we can apply equation (3. 1 ) to the shock front moving towards the right and evaluate 



it in the frame comoving with the contact discontinuity to obtain that the velocity ahead of the 
right propagating shock is 



V2 



{P3 -P2){e3' - 62) 
(e2 +P3){e3' +P2) 



(4.2) 



1 



Figure 2. Schematic wave pattern in the pressure for the decay of a discontinuity generating two shock 
waves propagating in opposite directions. The vertical solid lines show the position of the shock fronts 
while the dot-dashed vertical line shows the position of the contact discontinuity. The different arrows show 
the gas flow and the directions of propagation of the different fronts. 



Equations (4.1) and (1-.2) can now be used to derive the relativistic expression for the relative 
velocity of the flow ahead of the two shocks {vi2)2s- proved in Appendix A, the expression 
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for the relative velocity between the unperturbed states is a monotonic function of p3 for all 
possible wave patterns. In particular, for the present choice of initial data, this expression is a 
monotonically increasing function of p-^. As a result, the value of {vi2)^^ can be used to build a 
criterion for the occurrence of two shocks propagating in opposite directions. In fact, since pi is 
the smallest value that can take, two shocks will form if 



^ ^~ N _ {Pi -P2){e - 62) 



where 



hp ~ Pi = h- 



(e +P2)(e2 +Pi) ' 
IPi 



Pi , 



ii-m-i) 

and h is the only positive root of the Taub adiabat (Taub 1978, Maitf & Muller 1994) 



1 



(7 - 1)(P2 -P3) 



IPs. 



(7- 1)(P2 -Ps) r , h2{p2-pz) 



7P3 



P2 







(4.3) 



(4.4) 



(4.5) 



when p3 pi- Simple calculations reported in Appendix B show that the Newtonian limit of 
(■512)23 corresponds to the expression derived by Landau and Lifshitz (1987). 



5. LS^L^CR^TZ^R: One Shock and one Rarefaction Wave 

We next consider the wave pattern produced by one shock front propagating towards the right 
and one rarefaction wave propagating in the opposite direction (see Fig. ||). This situation is 
therefore characterized by pi > pa > P2. 
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Figure 3. Schematic wave pattern in the pressure for the decay of the discontinuity into a shock wave 
propagating towards the right and a rarefaction wave propagating in the opposite direction. The vertical 
lines show the discontinuities formed (continuous for the shock front; dashed for the head and the tail of 
the rarefaction wave; dot-dashed for the contact discontinuity) while the arrows show their direction of 
propagation and that of the gas flow. 



Evaluating expression (3.5) in the reference frame comoving with the contact discontinuity. 



we can evaluate the flow velocity ahead of the rarefaction wave to be 

1- A+(P3) 



Vl = 



l + A+{p^) 



(5.1) 
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where 



1)1/2 _c,(p3) 



(7-l)l/2+C,(p3) 



1)1/2 + c«(pi) 



(7-l)i/2-c.(pi) 



2/(7-1)1/^ 



(5.2) 



The flow velocity ahead of the shock front can be derived as in Section Q by evaluating equation 



( |3.1[ ) in the reference frame of the contact discontinuity to get 



W2 = 



' (P3 -P2)(e3' - 62) 

(es' +P2)(e2 +P3) ' 



(5.3) 



which, combined with expression (5.1 ), can be used to derive the relativistic expression for the 
relative velocity of the fluids ahead of the shock and ahead of the rarefaction wave (ui2)s^ . As 
for {vi2)^g, it can be shown that (wi2)s„ is a monotonically increasing function of (see Fig. |] 
and Appendix A for an analytic proof). Exploiting now the knowledge that for this wave pattern 
the pressure in the region between the two waves must satisfy p2 < Ps < Pi, we can establish 
that the criterion on the relative velocity for having one shock and one rarefaction wave is 



(^^12) 



l + A+{p3 



< V12 < 



' {Pi -P2){e - 62) 
(62 +pi)(e +P2) 



(^^12)2 



(5.4) 



Note that the upper limit of (5.4) coincides with (wi2)2S' which is the lower limit for the occur- 
rence of 2 shock waves (4.3) and whose Newtonian limit coincides with the equivalent one found 
by Landau and Lifshitz (1987) (see Appendix B). Note also that in the limit, — > p2, regions 1 
and 2 are connected by a single rarefaction wave. In this case the sound speed can be computed 
using p3 = p2 but with ~ pi{p2/p\Y^^ ■ Finally, note that this is the only wave pattern in 
which vi and V2 have the same sign and it therefore includes the classical shock-tube problem, 
where vi ^ V2 = 0. 



6. LTZ^L^CR-Jl^R: Two Rarefaction Waves 

We now consider the wave pattern produced by two rarefaction waves propagating in opposite 
directions (see Fig. This situation is characterized by pi > P2 > P3 and when the waves are 
sufficiently strong it might lead to a vacuum region (p^ = 0) behind the rarefaction waves. 





Figure 4. Schematic wave pattern in the pressure for the decay of the discontinuity into two rarefaction 
waves propagating in opposite directions. The vertical lines show the discontinuities formed (dashed for 
the head and the tail of the rarefaction waves; dot-dashed for the contact discontinuity) while the arrows 
show their direction of propagation and that of the gas flow. Note that the region downstream of the two 
rarefaction waves has a density ps > 0. 
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Following again the same procedure discussed in the previous Section, we can determine the 
values of the fluid velocities ahead of the two rarefaction waves as, respectively. 



Vi 



V2 = 



1-A^{py) 
1 + A-{P3') ' 



where 



A^{P3' 



(7-l)i/2 + c,(p3') 



(7-l)^/^ + c.(p2) 

1)1/2 



-2/(7-1)1/^ 



(6.1) 



(6.2) 



(6.3) 



_(^_ 1)1/2 „c,(p2) 

We have indicated with (ps' ) the sound speed in the region 3', which differs from the one in 
regi on 3 because of the jump in the densities pa and ps' . The relative velocity built using ( xl ) 
and ( 6^ ) is then 



(6.4) 



A+ips) + A^ipy) ■ 

As for the relative velocities of the previous wave patterns, it can be shown that (fi2)2H' 
a monotonically increasing function of p^ (see Fig. |l] and Appendix A for an analytic proof) so 
that the criterion for the occurrence of two rarefaction waves can be expressed as 



whe re A ^ {p^ 
for Q. 



P2) 



(^12)2fl|p3=o < ""12 < (Wl2)2n|p3^p^ " i^l2)sR , (6-5) 

= 1 and therefore the upper limit for (|6.5|) coincides with the lower limit 



(6.6) 



The condition (6.5) can also be expressed in a more useful form as 



^1-^2^ ^ l-A+{p2) 



Si + 5*2 

where, the constants Si and 5*2 are shorthand for 



l + A+{p2 



Si 



(7~l)^/^+c.(pi) 
(7-l)i/2-c.(pi) 

(7-l)^/^ + c,(p2) 
(7-l)i/2-c,(p2) 



2/(7-1)^/^ 



-2/(7-1)^/^ 



(6.7) 



(6.8) 



An important property of equation (6.4) is that it can be inverted analytically. This involves 
rewriting it in terms of a quartic equation in the unknown sound speed in the region (see 
Appendix C for the explicit form of the equation). Once the relevant real root of this equation has 
been calculated analytically, the value of the pressure can be found through a simple algebraic 
expression. In this case, therefore, the solution of the exact relativistic Riemann problem can be 
found in an analytic closed form. 

We conclude the analysis of the different wave patterns with a comment on the case of two 
rarefaction waves propagating in opposite directions and leaving behind them a region with zero 
density and pressure. This situation occurs when the fluids in the regions 1 and 2 are moving 
sufficiently fast in opposite directions. This is the case whenever the relative velocity between 
the two states ahead of the rarefaction waves is less than or equal to the lower hmit for (i'i2)2fj, 
i.e. 

Vl2 < {V12)2R = ^ ■ ^^'^^ 
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Also in this case, taking the Newtonian limit of (wi2)2b we obtain the corresponding expression 
derived by Landau and Lifshitz (1987) (see Appendix B). Finally, note that when a vacuum is 
produced, vi2 is no longer dependent on and this branch of the curve cannot be plotted in 




7. Numerical Implementation 

The core of most exact Riemann solvers, in both Newtonian and relativistic hydrodynamics, is 
based on the numerical computation of the pressure in the regions and i?* that form behind the 
waves. The key property exploited when performing the numerical calculation is that the velocity 
in such regions can be expressed as a monotonic function of the pressure i.e. u^, — u^, (pl,), 
and Ufl^ = vr^ (pr,)- Since there is no jump across a contact discontinuity in either the velocity 
or in the pressure, the numerical solution of the Riemann problem consists then of finding the 
root — PL, — PR, of the nonlinear equation 

vlAp*) - vrAp*) ^ , (7.1) 

where vl,,vr, have different functional forms according to the different wave patterns pro- 
duced. This method has two obvious disadvantages: (7) it cannot determine, using the initial 
conditions, the wave pattern produced and thus which of the functional forms to use for vl, , vr, ; 
f 2) it cannot provide a straightforward bracketing interval for the root. In practice, however, these 
difficulties are effectively balanced by efficient algorithms based on a sequence of trial and er- 
ror attempts that rapidly bracket the root and determine the correct equation to solve (see, for 
instance, Marti, & Miiller 1999 and the algorithm presented therein). 

The exact Riemann solver which we propose here differs from the one discussed above mostly 
because it avoids the disadvantages (]) and (2). In fact, as discussed in Section |l], by comparing 
the relative velocity between the initial left and right states (ui2)o with the relevant limiting val- 
ues constructed from the initial conditions (512)23' i^i2)sR' (^12)27?.' ^'^^ determine both the 
wave pattern which will be produced and the correct bracketing range in the pressure. Once this 
information has been obtained, the Riemann problem can be solved either through the solution of 



equation (7. 1 ) or, equivalently, by looking for the value of the pressure p» which would produce 
a relative velocity (wi2)o- This latter approach, which we will be discussing in the following, 
involves then the solution of the nonlinear equation 

wi2(p*) - ("12)0 = , (7.2) 

where Wi2(p*) is given by the expressions for (wi2)2S' or (^12)3/?' or (^12)2^ derived in Sec- 



tions ^|-|^. Furthermore, in the case of two rarefaction waves, equation (7.2) can also be solved 
analytically. 

Besides providing direct information about the wave pattern produced, about the correct equa- 
tion to solve and the relevant bracketing interval, our approach is also very simple to implement 
numerically. In practice, the basic steps for the solution of the Riemann problem can be sum- 
marised as follows: 

(fl) Evaluate from the initial conditions the three limiting relative velocities (■512)23' (^i2)sn' 

(b) Determine the wave pattern and the functional form of vi2{p*) by comparing (wi2)o with 
the Hmiting values calculated in (a) and according to the scheme below. 
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(c) According to the wave pattern found, determine the extremes Pmax and Pmin of the pres- 
sure interval bracketing . Within our conventions this is equivalent to setting^ 





s^cs^ 


Ln^LXR*S^R 


Ln^LXR*n^R 


pmin 


ma\{pi,p2) 


min(pi,p2) 





Pmax 


00 


max(pi,p2) 


min(pi,p2) 



id) Solve equation (7.2) and determine p*. 

(e) Complete the solution of the Riemann problem by computing the remaining variables of 
the intermediate states and i?* . 

We have implemented our algorithm for an exact Riemann solver and have tested it for a range 
of Riemann problems. We have also compared the performance of our algorithm with the "stan- 
dard" approach presented by Marti & Miiller 1999 and have found a systematic reduction in 
the computational costs for the same level of accuracy in the solution. The quantitative efficiency 
improvement depends on the type of problem under consideration. In the case of a generic hydro- 
dynamical problem (in which very simple Riemann problems are solved), our approach brackets 
the solution very closely and this produces substantial computational improvements of up to 
30%. However, for the specific cases discussed in this paper, where very strong shocks have 
been considered, the speed-up is smaller and of the order of 10%. It is worth noting that such an 
improvement could reduce appreciably the computational costs in three-dimensional relativistic 
hydrodynamics codes, where this operation is repeated a very large number of times. 



8. Conclusion 

We have presented a new procedure for the numerical solution of the exact Riemann problem 
in relativistic hydrodynamics. In this approach special attention is paid to the relativistic invariant 
expression for the relative velocity V12 between the unperturbed left and right states. This has 
been shown to be a monotonic function of the pressure in the region formed between the wave 
fronts. The determination of this pressure is the basic step in the solution of the exact Riemann 
problem. 

The use of the relative velocity has a number of advantages over alternative exact Riemann 
solvers discussed in the literature. In particular, it extracts the information implicitly contained 
in the data for two initial states to deduce the wave pattern that will be produced by the decay of 
the discontinuity between these two states. This, in turn, allows an "a priori" determination to be 
made of the interval in pressure bracketing and the correct functional form for the nonlinear 

I Note that in practice, the upper hmit for the pressure in the case of two shocks is found by starting 
from a reasonable value above Pmin, which is incremented until the solution is effectively bracketed. 
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equation whose root will solve the exact Riemann problem. All of these advantages translate, 
in practice, into a simpler algorithm to implement and an improved efficiency in the numerical 
solution of the Riemann problem. Furthermore, in the case of two rarefaction waves propagating 
in opposite directions (and in strict analogy with what happens in Newtonian hydrodynamics), 
the use of the relative velocity allows for the solution of the Riemann problem in an analytic 
closed form. 

Because of all of the advantages discussed above, its intrinsic simplicity of implementation and 
the numerical efficiency gain it produces, this new exact Riemann solver should be considered 
as an interesting alternative to the traditional exact Riemann solver presently discussed in the 
literature. Investigations about the extension of this approach to multidimensions are currently in 
progress. 

It is a pleasure to thank J. A. Font, C. Gheller, and J. C. Miller for useful discussions. Finan- 
cial support for this research has been provided by the Italian Ministero dell' Universita e della 
Ricerca Scientifica and by the EU Programme 'Improving the Human Research Potential and 
the Socio-Economic Knowledge Base' (Research Training Network Contract HPRN-CT-2000- 
00137)." 



Appendix A: Monotonicity of the relative velocity as function of 

We here prove that vi2 is a monotonic function of for all of the possible wave patterns. In 
particular, because of our choice of referring the initial left state as to "1" and the initial right 
state as to "2", we will here show that vi2 is a monotonically increasing function of p*. 

Denoting by a' the first derivative of the quantity a with respect to p* = = py, it is 
straightforward to obtain that the first derivative of ( 2.5 ) is 



V[{l-Vl)-V'2{l-Vl) 
(1 - V1V2Y 



(8.1) 



Since vi < 1, V2 < 1, the two terms in the parentheses of (B.l) are always positive, and the 
proof that V12 is monotonically increasing will follow if it can be shown that v'-^ and — are both 
positive. We wiU do so for each of the three possible wave patterns. 



Two Shock Fronts 

Taking the derivative of the fluid velocity ahead of the left propagating shock front (measured 



from the contact discontinuity) [cf. equation ( 4-. 1 )] we obtain 



(ei [(63 - ei)(e3 + (p3 -pi)(ei +P3)e3] 

2wi(ei +P3)^(e3 +PiY 



(8.2) 



Since the energy density is an increasing function of pressure, 63 > 0; furthermore, p-^ > pi and 
63 > ei for the wave pattern considered and v[ > as a result. The equivalent expression for 
the derivative of the fluid velocity ahead of the right propagating shock front ( meas ured from the 
contact discontinuity) [cf. eq (4-.2)] can be obtained by replacing in equation (8.2) the indices 1 
and 3 by 2 and 3' respectively, i.e. 



(62 +P2) [(63' - e2){e3> +P2) + {P3' - P2){e2 + P3')e3,] 
2^2(62 +P3')^(63' +P2)^ 



(8.3) 



Since now for the wave pattern considered: V2 < 0, ps' > p2 and 63/ > 62, we are led to 
conclude that —V2 > 0, thus making the overall v'12 positive for any value of p^. 
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One Shock and one Rarefaction Wave 

In this case we only need to show that v'i> Q since for the velocity ahead of the right propagating 
shock front we can use the results derived in (8.3). Taking the derivative of expression (5.1 ) then 
yields 

2a;(p3) 



[l + ^+(P3)]^ 



(8.4) 



where 



1)1/2 



(^-1)1/2 
^ -C14(P3) , 



Cs{pi) 



^+(P3) 



(2-v^)/2 



[(^_ 1)1/2 +e.(p3)]^ 



(8.5) 



where Ci > 0. When the sound speed Cs (pa) is an increasing function of pressure, as is the case 
for a polytropic equation of state [cf. eq (|2.2[)], v[> Q and therefore v'12 > 0. 



Two Rarefaction Waves 

What we need to show in this case is that ^2 < since we can exploit the previous result that 
v'i> {) where vi is the fluid velocity ahead of the left propagating rarefaction wave as measured 
from the contact discontinuity. In this case, taking the derivative of expression (6.2) one obtains 



(8.6) 



where now 



A'_(p3') = 4 



(7-l)l/2-c.(p2) 



A_(p3,)(2+y^)/2 



[(7-l)V2+C,(p3')]- 



(8.7) 



with C2 > and therefore V2 < 0. 

We have therefore shown that for all of the wave patterns considered v[ > and —V2 > 0, 
thus proving that V12 is always a monotonically increasing function of p,. 



Appendix B: Newtonian limits of {1)12)23^ ("^12)3^? (^12)2^ 

We here show that the three limiting values of (■012)25' (''^i2)sh ^^'^ (^12)2/! reduce to their 
Newtonian counterparts in the limit of w, Cs 0, and h ^ 1. In particular, we will restrict 



ourselves to considering the case of a polytropic equation of state (2.2). 

We start by considering the Newtonian limit of (wi2)2s which is obtained when p/e ^ 1 and 
e l/V, with V ^ 1/p being the specific volume. In this case, then 



(^'12)2 



Nowt 



{pi - P2)(e " 62) 

662 



= Vbl-P2) (1/62-1/6). (8.1) 



which coincides with the corresponding expression derived by Landau and Lifshitz (1987) but 
with inverted indices. 

We next consider the Newtonian limit of {V12) sr which is obtained when both 65(^3) ^ 1 
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and Cs(pi) ^ 1. In this case, 

A+{p^) ^ ( 1 



1 



Cs(p3) 

4c, (p3) 
7-1 



Cs(pi) 



7-1 



7-1 



(Cs(p3) - Cs[pi)) 



(8.2) 



so that the Newtonian limit of (wi2)sn is given by 

_ 1-A+(P3) 



(Wl2)c 



Newt 



7-1 



Cs(Pl) 



Cs{pi) 



Bearing in mind that 



Cs(Pl) 



(7-l)/27 



we obtain 



(«12)sH 



-Cs(Pl) 



1 



(7-l)/27' 



(8.3) 



(8.4) 



(8.5) 



Newt T 

which again coincides with the corresponding expression derived by Landau and Lifshitz (1987) 
but with inverted indices. 

Finally, we consider the Newtonian limit of ("512 ) 2 j, for Cs(pi), Cs(p2) <C 1. In this case 



S2 ^ 1 



4c, (pi) 
7-1 

4c, (P2) 
7-1 



(8.6) 



so that the Newtonian limit is given by 



(^12) 2R 



_ S2 — Si 

Newt + 'S'2 



4c,(pi) _ 4c,(p2) 
7 — 1 7 — 1 

2cs(pi) _ 2c, (P2) 
7 — 1 7 — 1 
2c, (pi) _ 2c, (pa) 
7 — 1 7 — 1 



n -1 



7-1 
2 

7^ 



[c,(pi) - Cs(p2)] 
[C,(pi) - C,(p2)] 



(8.7) 



Once more, expression (8.7) coincides with the corresponding expression derived by Landau and 
Lifshitz (1987) but with inverted indices. 



Appendix C: A closed form solution in the case of Two Rarefaction Waves 

As discussed in Sections^ and ^ when (wi2)2H < (^12)0 < (i3i2)sn, the initial conditions give 
rise to two rarefaction waves and it is possible to derive a closed form solution for the unknown 
pressure p* . In this way we can, at least in principle, avoid any numerical root finding procedure 
and determine the solution exactly. In this Appendix we first derive this analytic solution in the 
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context of relativistic hydrodynamics and then calculate its Newtonian limit. We will restrict 
ourselves to considering the particular case of a polytropic equation of state (2.2). 

Using expression (3.4) we can write the pressures and p^' as functions of the sound speeds 
Cs {p3 ) and Cs {ps' ) which, for convenience, we will hereafter refer to as x and x' respectively 

17/(7-1) 



P3 



P-3' 



-1/(7-1) [ ^ (7-1) 

7(7 — 1) — 72:^ 



7/(7-1) 



(8.1) 



(8.2) 



.7(7 - 1) - 7(2:')^ 

where fci = pi/pj and fc2 = P2I pI the two constants entering the polytropic equation. Since 
Vs = Ps' = P*, we can obtain the following relation between x' and x: 

/^2_ 7(7-1)2^^ 



7x^(1 — a) + a7(7 — 1) 



(8.3) 



where a = (fci/fc2)^^^- The expression for the relative velocity (6.4) can also be written as 

A+jps) _ 1 - {vi2)o 
A-{p3') 1 + (^12)0 



(8.4) 



and we then use expressions ( p.2[ ) and ( |6.3| ) to expand the left hand side of ( |8.4| ). After some 
algebra we are then left with 

("12)0 



r - Cs{p2) 



r + Cs{p2) 



r - Cs{pi) 



r + Cs(pi) 



1 



1 + (^12)0 



where = 7 — 1 and the right hand side of (8.5) is a constant which we rename as 



n = 



" r - Cs(p2) 

T + Cs{p2) 

Introducing now the auxiliary quantity 



r - Csipi) 



r + Cs{pi) 



1 - (1^12)0 



1 + (l'l2)o 



r/2 



/3 = 



1 + n 
1 -n 



expression (8.5) can be written as 



l\2 



{x') 



T{xp - r) 

x-Tf3 



(8.5) 



(8.6) 



(8.7) 



(8.8) 



Comparing now (3.3 ) and (8.8) gives a 4-th order equation in the unknown sound velocity x 



qqx + aix + a2X + a^x + 04 = , 



where 



ao = 1-/3^(1 -a) 
a2=r\l-a){(3^ - 



1) 



04 



(8.9) 

(8.10) 
(8.11) 
(8.12) 
(8.13) 
(8.14) 



The analytic solution of equation (8.9) will yield at least two real roots, one of which will be the 
physically acceptable one: i.e. positive, less than one, and such that the pressure falls in the 
relevant bracketing interval. 
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In its Newtonian limit, equation ( |8.9| ) reduces to a second order equation in the unknown sound 
velocity 

^^-l^a;2 + 2i;x-i;2 = 0, (8.15) 

where 

E = c,(pi) + c,(p2) + -^^^^y^i^i2 , (8.16) 

and vi2 = vi — V2- The fact that the Newtonian Riemann problem in the case of two rarefac- 
tion waves can be solved analytically is well known and is at the basis of the so called "Two 
Rarefaction Approximate Riemann Solver" (Toro 1997). 



REFERENCES 

COLELLA, P. AND WOODWARD, P. R. 1984, The Piecewise Parabolic Method (PPM) for Gas-Dynamical 

Simulations, / Comput. Phys., 54, 174. 
COURANT, R. AND Friedrichs, K. O. 1948, Supersonic Flows and Shock Waves, Interscience. 
Gheller, C. 1997, A High Resolution Hydrodynamical Code for the Study of Cosmological Structures, 

PhD Thesis, SISSA. 

GODUNOV, S. K. 1959, A Finite Difference Method for the Numerical Computation and Discontinuous 

Solutions of the Equations of Fluid Dynamics, Mat. Sb., 47, 271. 
IBANEZ, J.M. & Mart!, J.M. 1999, Riemann Solvers in Relativistic Astrophysics, / Comput. Appl. 

Math., in press. 

Landau, L. D. and Lifshitz, E. M. 1987, Fluid Mechanics (Second Edition), Pergamon. 

Marti', J.M. & Muller, E. 1994, The Analytical Solution of the Riemann Problem in Relativistic 

Hydrodynamics, / Fluid Meek., 258, 317. 
Mart!, J.M. & MOller, E. 1996, Extension of the Piecewise Parabohc Method to One-Dimensional 

Relativistic Hydrodynamics, / Comput. Phys. 123, 1. 
Mart!, J.M. & MOll er, E. 1999 Numerical Hydrodynamics in Special Relativity, Living Reviews in 



General Relativity, astro-ph/ 9906333 



Pons, J. A., Font, J. A., Marti, J.M., Ibanez, J.M. & Miralles, J. A. 1998, General Relativistic 
Hydrodynamics with Special Relativistic Riemann Solvers, Astron. Astroph. 339, 638. 

Pons, J. A., Marti, J.M., MOller, E. 2000, The exact solution of the Riemann problem with non zero 
tangential velocities in relativistic hydrodynamics J. Fluid Mech. 422, 125. 

Roe, p. L. 1981, Approximate Riemann Solvers, Parameter Vectors and Difference Scheme, J. Comp. 
Phys., 43, 357. 

Sod, G. a. 1978, A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic 

Conservation Laws, /. Comp. Phys., 27, 1. 
Taub, a. H., 1948, Relativistic Rankine-Hugoniot Relations, Phys. Rev, 74, 328. 
Taub, a. H., 1978, Relativistic Fluid Mechanics, Ann. Rev Fluid Mech., 10, 301. 
Toro, E. 1991 , Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer. 



